//  Copyright (c) 2011 M.A. (Thijs) van den Berg, http://sitmo.com/
//
//  Use, modification and distribution are subject to the MIT Software License. 
// (See accompanying file LICENSE.txt or copy at http://www.bmlib.org/LICENSE.txt)


#ifndef BMLIB_GBM_OPTION_BARRIER_HPP
#define BMLIB_GBM_OPTION_BARRIER_HPP

#include <cmath>
#include "math_normal.hpp"

namespace bmlib{

namespace gbm_option_barrier_detail {
    
    enum BarrierTypes
    {
        barriertype_DIC,
        barriertype_DIP,
        barriertype_DOC,
        barriertype_DOP,
        barriertype_UIC,
        barriertype_UIP,
        barriertype_UOC,
        barriertype_UOP
    };
    
    inline double price(
            double S, double Y, double v,   // underlying gbp
            double r,                       // rate for NPV
            double X, double H, double t,   // option contract parameters
            BarrierTypes BarrierType
            )
    {
    
        // common terms
        double K=1.;
        double v_sqrt_t = v*sqrt(t);
        
        double u = (Y-v*v/2.0)/v/v;
        double lambda = sqrt(u*u + 2.0*r/v/v);
        double x1 = log(S/X)    /v_sqrt_t + (1.0+u)*v_sqrt_t;
        double x2 = log(S/H)    /v_sqrt_t + (1.0+u)*v_sqrt_t;
        double y1 = log(H*H/S/X)/v_sqrt_t + (1.0+u)*v_sqrt_t;
        double y2 = log(H/S)    /v_sqrt_t + (1.0+u)*v_sqrt_t;
        double z  = log(H/S)    /v_sqrt_t +  lambda*v_sqrt_t;
        double nu,o;
        
        switch (BarrierType) {
            case barriertype_DIC:
            {
                nu = 1; o = 1;
                double A = o*S*exp((Y-r)*t)*N(o*x1) - o*X*exp(-r*t)*N(o*x1-o*v*sqrt(t));
                double B = o*S*exp((Y-r)*t)*N(o*x2) - o*X*exp(-r*t)*N(o*x2-o*v*sqrt(t));
                double C = o*S*exp((Y-r)*t)*pow(H/S,2.0*(u+1.0))*N(nu*y1) - o*X*exp(-r*t)*pow(H/S,2.0*u)*N(nu*y1-nu*v*sqrt(t));
                double D = o*S*exp((Y-r)*t)*pow(H/S,2.0*(u+1.0))*N(nu*y2) - o*X*exp(-r*t)*pow(H/S,2.0*u)*N(nu*y2-nu*v*sqrt(t));
                double E = K*exp(-r*t)*( N(nu*x2-nu*v*sqrt(t)) - pow(H/S,2.0*u)*N(nu*y2-nu*v*sqrt(t)));
                double F = K*( pow(H/S,u+lambda)*N(nu*z) + pow(H/S,u-lambda)*N(nu*z-2.0*nu*lambda*v*sqrt(t)) );                
                if (X>H)    return C + E; 
                else        return A - B + D + E;
                break;
            }
            case barriertype_UIC:
            {
                nu = -1; o = 1;
                double A = o*S*exp((Y-r)*t)*N(o*x1) - o*X*exp(-r*t)*N(o*x1-o*v*sqrt(t));
                double B = o*S*exp((Y-r)*t)*N(o*x2) - o*X*exp(-r*t)*N(o*x2-o*v*sqrt(t));
                double C = o*S*exp((Y-r)*t)*pow(H/S,2.0*(u+1.0))*N(nu*y1) - o*X*exp(-r*t)*pow(H/S,2.0*u)*N(nu*y1-nu*v*sqrt(t));
                double D = o*S*exp((Y-r)*t)*pow(H/S,2.0*(u+1.0))*N(nu*y2) - o*X*exp(-r*t)*pow(H/S,2.0*u)*N(nu*y2-nu*v*sqrt(t));
                double E = K*exp(-r*t)*( N(nu*x2-nu*v*sqrt(t)) - pow(H/S,2.0*u)*N(nu*y2-nu*v*sqrt(t)));
                double F = K*( pow(H/S,u+lambda)*N(nu*z) + pow(H/S,u-lambda)*N(nu*z-2.0*nu*lambda*v*sqrt(t)) );                
                if (X>H)    return A + E; 
                else        return B - C + D + E;
                break;
            }
            case barriertype_DIP:
            {
                nu = 1; o = -1;
                double A = o*S*exp((Y-r)*t)*N(o*x1) - o*X*exp(-r*t)*N(o*x1-o*v*sqrt(t));
                double B = o*S*exp((Y-r)*t)*N(o*x2) - o*X*exp(-r*t)*N(o*x2-o*v*sqrt(t));
                double C = o*S*exp((Y-r)*t)*pow(H/S,2.0*(u+1.0))*N(nu*y1) - o*X*exp(-r*t)*pow(H/S,2.0*u)*N(nu*y1-nu*v*sqrt(t));
                double D = o*S*exp((Y-r)*t)*pow(H/S,2.0*(u+1.0))*N(nu*y2) - o*X*exp(-r*t)*pow(H/S,2.0*u)*N(nu*y2-nu*v*sqrt(t));
                double E = K*exp(-r*t)*( N(nu*x2-nu*v*sqrt(t)) - pow(H/S,2.0*u)*N(nu*y2-nu*v*sqrt(t)));
                double F = K*( pow(H/S,u+lambda)*N(nu*z) + pow(H/S,u-lambda)*N(nu*z-2.0*nu*lambda*v*sqrt(t)) );                
                if (X>H)    return B - C + D + E;
                else        return A + E; 
                break;
            }
            case barriertype_UIP:
            {
                nu = -1; o = -1;
                double A = o*S*exp((Y-r)*t)*N(o*x1) - o*X*exp(-r*t)*N(o*x1-o*v*sqrt(t));
                double B = o*S*exp((Y-r)*t)*N(o*x2) - o*X*exp(-r*t)*N(o*x2-o*v*sqrt(t));
                double C = o*S*exp((Y-r)*t)*pow(H/S,2.0*(u+1.0))*N(nu*y1) - o*X*exp(-r*t)*pow(H/S,2.0*u)*N(nu*y1-nu*v*sqrt(t));
                double D = o*S*exp((Y-r)*t)*pow(H/S,2.0*(u+1.0))*N(nu*y2) - o*X*exp(-r*t)*pow(H/S,2.0*u)*N(nu*y2-nu*v*sqrt(t));
                double E = K*exp(-r*t)*( N(nu*x2-nu*v*sqrt(t)) - pow(H/S,2.0*u)*N(nu*y2-nu*v*sqrt(t)));
                double F = K*( pow(H/S,u+lambda)*N(nu*z) + pow(H/S,u-lambda)*N(nu*z-2.0*nu*lambda*v*sqrt(t)) );                
                if (X>H)    return A - B + D + E;
                else        return C + E; 
                break;
            }
            case barriertype_DOC:
            {
                nu = 1; o = 1;
                double A = o*S*exp((Y-r)*t)*N(o*x1) - o*X*exp(-r*t)*N(o*x1-o*v*sqrt(t));
                double B = o*S*exp((Y-r)*t)*N(o*x2) - o*X*exp(-r*t)*N(o*x2-o*v*sqrt(t));
                double C = o*S*exp((Y-r)*t)*pow(H/S,2.0*(u+1.0))*N(nu*y1) - o*X*exp(-r*t)*pow(H/S,2.0*u)*N(nu*y1-nu*v*sqrt(t));
                double D = o*S*exp((Y-r)*t)*pow(H/S,2.0*(u+1.0))*N(nu*y2) - o*X*exp(-r*t)*pow(H/S,2.0*u)*N(nu*y2-nu*v*sqrt(t));
                double E = K*exp(-r*t)*( N(nu*x2-nu*v*sqrt(t)) - pow(H/S,2.0*u)*N(nu*y2-nu*v*sqrt(t)));
                double F = K*( pow(H/S,u+lambda)*N(nu*z) + pow(H/S,u-lambda)*N(nu*z-2.0*nu*lambda*v*sqrt(t)) );                
                if (X>H)    return A - C +  F;
                else        return B - D +  F; 
                break;
            }
            case barriertype_UOC:
            {
                nu = -1; o = 1;
                double A = o*S*exp((Y-r)*t)*N(o*x1) - o*X*exp(-r*t)*N(o*x1-o*v*sqrt(t));
                double B = o*S*exp((Y-r)*t)*N(o*x2) - o*X*exp(-r*t)*N(o*x2-o*v*sqrt(t));
                double C = o*S*exp((Y-r)*t)*pow(H/S,2.0*(u+1.0))*N(nu*y1) - o*X*exp(-r*t)*pow(H/S,2.0*u)*N(nu*y1-nu*v*sqrt(t));
                double D = o*S*exp((Y-r)*t)*pow(H/S,2.0*(u+1.0))*N(nu*y2) - o*X*exp(-r*t)*pow(H/S,2.0*u)*N(nu*y2-nu*v*sqrt(t));
                double E = K*exp(-r*t)*( N(nu*x2-nu*v*sqrt(t)) - pow(H/S,2.0*u)*N(nu*y2-nu*v*sqrt(t)));
                double F = K*( pow(H/S,u+lambda)*N(nu*z) + pow(H/S,u-lambda)*N(nu*z-2.0*nu*lambda*v*sqrt(t)) );                
                if (X>H)    return F;
                else        return A - B + C - D + F; 
                break;
            }
            case barriertype_DOP:
            {
                nu = 1; o = -1;
                double A = o*S*exp((Y-r)*t)*N(o*x1) - o*X*exp(-r*t)*N(o*x1-o*v*sqrt(t));
                double B = o*S*exp((Y-r)*t)*N(o*x2) - o*X*exp(-r*t)*N(o*x2-o*v*sqrt(t));
                double C = o*S*exp((Y-r)*t)*pow(H/S,2.0*(u+1.0))*N(nu*y1) - o*X*exp(-r*t)*pow(H/S,2.0*u)*N(nu*y1-nu*v*sqrt(t));
                double D = o*S*exp((Y-r)*t)*pow(H/S,2.0*(u+1.0))*N(nu*y2) - o*X*exp(-r*t)*pow(H/S,2.0*u)*N(nu*y2-nu*v*sqrt(t));
                double E = K*exp(-r*t)*( N(nu*x2-nu*v*sqrt(t)) - pow(H/S,2.0*u)*N(nu*y2-nu*v*sqrt(t)));
                double F = K*( pow(H/S,u+lambda)*N(nu*z) + pow(H/S,u-lambda)*N(nu*z-2.0*nu*lambda*v*sqrt(t)) );                
                if (X>H)    return A - B + C - D + F;
                else        return F; 
                break;
            }
            case barriertype_UOP:
            {
                nu = -1; o = -1;
                double A = o*S*exp((Y-r)*t)*N(o*x1) - o*X*exp(-r*t)*N(o*x1-o*v*sqrt(t));
                double B = o*S*exp((Y-r)*t)*N(o*x2) - o*X*exp(-r*t)*N(o*x2-o*v*sqrt(t));
                double C = o*S*exp((Y-r)*t)*pow(H/S,2.0*(u+1.0))*N(nu*y1) - o*X*exp(-r*t)*pow(H/S,2.0*u)*N(nu*y1-nu*v*sqrt(t));
                double D = o*S*exp((Y-r)*t)*pow(H/S,2.0*(u+1.0))*N(nu*y2) - o*X*exp(-r*t)*pow(H/S,2.0*u)*N(nu*y2-nu*v*sqrt(t));
                double E = K*exp(-r*t)*( N(nu*x2-nu*v*sqrt(t)) - pow(H/S,2.0*u)*N(nu*y2-nu*v*sqrt(t)));
                double F = K*( pow(H/S,u+lambda)*N(nu*z) + pow(H/S,u-lambda)*N(nu*z-2.0*nu*lambda*v*sqrt(t)) );                
                if (X>H)    return B - D + F;
                else        return A - C + F; 
                break;
            }
        }
    }
    
} // namespace gbm_option_barrier_detail



inline double gbm_barrier_down_and_in_call_price(double S, double Y, double v, double r, double X , double H, double t)
    { return gbm_option_barrier_detail::price(S, Y, v, r, X, H, t, gbm_option_barrier_detail::barriertype_DIC); }

inline double gbm_barrier_down_and_in_put_price(double S, double Y, double v, double r, double X , double H, double t)
    { return gbm_option_barrier_detail::price(S, Y, v, r, X, H, t, gbm_option_barrier_detail::barriertype_DIP); }

inline double gbm_barrier_down_and_out_call_price(double S, double Y, double v, double r, double X , double H, double t)
    { return gbm_option_barrier_detail::price(S, Y, v, r, X, H, t, gbm_option_barrier_detail::barriertype_DOC); }

inline double gbm_barrier_down_and_out_put_price(double S, double Y, double v, double r, double X , double H, double t)
    { return gbm_option_barrier_detail::price(S, Y, v, r, X, H, t, gbm_option_barrier_detail::barriertype_DOP); }

inline double gbm_barrier_up_and_in_call_price(double S, double Y, double v, double r, double X , double H, double t)
    { return gbm_option_barrier_detail::price(S, Y, v, r, X, H, t, gbm_option_barrier_detail::barriertype_UIC); }

inline double gbm_barrier_up_and_in_put_price(double S, double Y, double v, double r, double X , double H, double t)
    { return gbm_option_barrier_detail::price(S, Y, v, r, X, H, t, gbm_option_barrier_detail::barriertype_UIP); }

inline double gbm_barrier_up_and_out_call_price(double S, double Y, double v, double r, double X , double H, double t)
    { return gbm_option_barrier_detail::price(S, Y, v, r, X, H, t, gbm_option_barrier_detail::barriertype_UOC); }

inline double gbm_barrier_up_and_out_put_price(double S, double Y, double v, double r, double X , double H, double t)
    { return gbm_option_barrier_detail::price(S, Y, v, r, X, H, t, gbm_option_barrier_detail::barriertype_UOP); }


}
#endif